tic

%%% This is the mail file to solve the discrete time model.

clear all
close all

load initialV.mat

global parms i Sf kf Nf dVk_f dVS T0

y0 = 1.0;     % output of goods y0 - choose units so this is 1.0
%
Sbar = 2126.0527;    % Total feasible resources
Res_0 = 1./0.065100642;       % Initial estimate of proved reserves
ratio_0 = Sbar - Res_0;    % Initial ratio of Alpha2 to Alpha3
g_0 = 0.11068716;      % Initial value of per unit mining cost
n0 = 0.00832913;      % Initial level of mining investment
%
% % gN_0 = -0.081639946;    % Partial derivative of g with respect to N at t=0
% % Alpha3 = -gS_0.*Rho_0./gN_0;
%
alpha3 = 15; %-gS_0.*Rho_0./gN_0;
alpha2 = ratio_0.*alpha3;
%
delta = 0.04;
%
c0 = 0.6619974; % Initial value of consumption for calculate lambda_0
k0 = 3.6071282734; % Initial value of K for the differential equation
%
A =  y0./k0;
%
% Marginal cost of backstop energy p = (Gamma1+H)^(-alpha)
% Alpha is the slope of the learning curve
%
alpha = 0.25;
%
% Given alpha, Gamma1 determines the initial cost of renewable energy. Here we set it
% to 4 times the initial cost of fossil fuel.
%
Gamma1 = (4*g_0).^(-1./alpha);
%
% % After some time t, marginal cost will decline to Gamma2 and remain there.
% % We assume this ultimate minimum marginal cost of renewables
% % is 20% of the initial cost p when H = 0:
%
Gamma2 =0.8*g_0;
%
Abar = A*(1-Gamma2)+ (1-delta);
% Note: We should have Abar > 0. It is different from Abar in continuous
% time model

% To make discretized model comparable to continuous one, beta is redefined
% so that the long term growth rate of the two models are the same.
% (beta*Abar)^(1/gamma)-1 = 4.07%, from which beta = 0.9673
beta  =  0.9673;

% psi is the effect of learning by doing on H
% psi has to be between 0 and 1
% A smaller value of psi will allow a larger role for explicit investment
% in renewable technology as opposed to learning by doing

psi = 0.33;
psr = 1/psi;
psc = 1-psi;
%
% Q = population growth that is used to "scale" resource extraction
% (variables other than fossil fuel exploitation are in per capita terms)
% popgr is the exogenous population growth rate

Q0 = 1;
popgr = 0.01;
%
% % gamma is the coefficient of relative rsik aversion. If we wish
% % to calibrate to a particular initial consumption level, we can
% % allow gamma to vary.
gamma = 4.0;
gamc = 1-gamma;
gamr = 1/gamma;

gS_0 = 0.00015;  % Partial derivative of g with respect to S at t=0
alpha1 = gS_0.*Res_0.^2;
alpha0 = g_0 - alpha1./Res_0;

parms = [delta A Sbar alpha0 alpha1 alpha2 alpha3 Gamma1 alpha psi beta ...
    gamma Q0 popgr Gamma2 Abar];
% % Note that parms is a 1x14 vector



%%%%%%%%%%%%%%%  Regime 1: B>0, j>=0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% A very preliminary algorithm:
%%% 1. Given Tbar, decrease S_T0 first, until S_T0 is too small to solve the
%%% problem appropriately. Note that S_T0 has little effects on k0, but
%%% decreases N0 and S0, especially S0.
%%% 2. Increase k_Tbar until it is too large to solve the problem, which will
%%% increase k0, while decrease both N0 and S0.
%%% 3. If k0 is still smaller than the target, while N0 and S0 positive,
%%% decrease Tbar and repeat step 1 and 2 again.


%%% solution for guesses in continuous model
% Tbar = 315.8;
% k_Tbar = 141704.98998437249;
% S_T0 = 1613;
%%%
%Guess the time Tbar when the economy transits to the analytical
% regime. Tbar has no effects on renewable regime.
 % integer
% % Tbar = 315;

% 2. Guess the capital stock at k_Tbar. k_Tbar cannot be smaller than some
% value between [140938.6464,140938.6465], or H can never hit zero. On the
% other hand, the larger the k_Tbar is, the shorter the renewable regime
% is. k_T0 will be larger as well.
% k_Tbar has to large than .648948 to make T0 no less than 88.
% % k_Tbar = 140938.648945;

% 3. Guess S_T0 = 1550; Sbar = 2126.0527
Tbar = 304;
k_Tbar = 140938.6510359;
S_T0 = 1582.22;
% % S_T0 = 1504;


% ======================================
% At Tbar, the marginal cost of backstop technology becomes Gamma2. Hence
% we have (Gamma1+H)^(-alpha)=Gamma2
H_Tbar = Gamma2^(-1/alpha)-Gamma1;

% Calculate k_t, H_t and V_t for all t in the renewable regime.

% Initialize the solution matrix
tt = 300; % tt should be larger than # of years in renewable regime
Vr = zeros(tt,1);
ir = zeros(tt,1);
jr = zeros(tt,1);
cr = zeros(tt,1);
kr = zeros(tt,1);
Hr = zeros(tt,1);
dVk = zeros(tt,1);
dVH = zeros(tt,1);

% Value of variables at Tbar
kr(1) = k_Tbar;
Hr(1) = H_Tbar;
ir(1) = ((beta*Abar)^gamr-1+delta)*k_Tbar;
cr(1) = (1-Gamma2)*A*kr(1)-ir(1);
dVk(1) = Abar*cr(1)^(-gamma);
dVH(1) = cr(1)^(-gamma)*(alpha*A*kr(1)*(Gamma1+Hr(1))^(-alpha-1));
i = 2;
while true
    
    
    % When i = 2, we are solving period Tbar-1. H = H_Tbar and j=0 at Tbar+1, djk and
    % djH at Tbar = 0.
    % Htilda = Hr(i-1); does not affect the results
    
    Cbar = psc*dVH(i-1)/dVk(i-1);
    cr(i) = (beta*dVk(i-1)).^(-gamr);
    
    j = @(k) Cbar.^psr*A*k;
    inv = @(k) kr(i-1)-(1-delta)*k;
    H = @(k) Hr(i-1)-Cbar^(psr*psc)*A*k;
    
    if H(kr(i-1))>0 && inv(kr(i-1))>0
        kr(i) = fzero(@(k) (1-(Gamma1+H(k)).^(-alpha))*A.*k-inv(k)-j(k)-cr(i),kr(i-1));
        Hr(i) = H(kr(i));
        
        if Hr(i) == Hr(i-1)
            error('H converges to a positive number')
        end
        jr(i) = j(kr(i));
        ir(i) = inv(kr(i));
        Vr(i) = cr(i).^gamc/gamc+beta*Vr(i-1);
        dVk(i) = cr(i)^(-gamma)*(A-(Gamma1+Hr(i))^(-alpha)*A)+...
            beta*(1-delta)*dVk(i-1)+beta*psi*dVH(i-1)*A^psi*kr(i)^(-psc)*jr(i)^psc;
        dVH(i) =  cr(i)^(-gamma)*(alpha*A*kr(i)*(Gamma1+Hr(i))^(-alpha-1))+beta*dVH(i-1);
        i = i+1;
    else
        
        Hr(i) = 0;
        kr(i) = fzero(@(k) (1-Gamma1.^(-alpha))*A.*k-inv(k)-j(k)-cr(i),kr(i-1));
        jr(i) = j(kr(i));
        ir(i) = inv(kr(i));
        Vr(i) = cr(i).^gamc/gamc+beta*Vr(i-1);
        dVk(i) = cr(i)^(-gamma)*(A-(Gamma1+Hr(i))^(-alpha)*A)+...
            beta*(1-delta)*dVk(i-1)+beta*psi*dVH(i-1)*A^psi*kr(i)^(-psc)*jr(i)^psc;
        dVH(i) =  cr(i)^(-gamma)*(alpha*A*kr(i)*(Gamma1+Hr(i))^(-alpha-1))+beta*dVH(i-1);
        % Cut out the zero lines of solution matrix
        yearsr = i;
        Hr = Hr(1:yearsr);
        kr = kr(1:yearsr);
        jr = jr(1:yearsr);
        ir = ir(1:yearsr);
        Vr = Vr(1:yearsr);
        cr = cr(1:yearsr);
        
           break
    end
    
end

T0 = Tbar - yearsr;

% Compare the results of the two models.
figure(1)
subplot(2,2,1), plot(T0+1:Tbar,kr(end:-1:1))
hold on
plot (T0_0+1:Tbar_0,kr0,'red')
title('k')

subplot(2,2,2),  plot(T0+1:Tbar,Hr(end:-1:1))
hold on
plot (T0_0+1:Tbar_0,Hr0,'red')
title('H')

subplot(2,2,3),  plot(T0+1:Tbar,ir(end:-1:1))
hold on
plot (T0_0+1:Tbar_0,ir0,'red')
title('i')

subplot(2,2,4),  plot(T0+1:Tbar,jr(end:-1:1))
hold on
plot (T0_0+1:Tbar_0,jr0,'red')
title('j')


%%%%%%%%%%%%%%%  Regime 2: R>0, n>=0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Key state variables: k and N
% Control variables: i and n
% Other state variables: Q and S
% Transition point T0

% Set the initial results
Vf = zeros(T0+1,1);
invf = zeros(T0+1,1);
nf = zeros(T0+1,1);
kf = zeros(T0+1,1);
Nf = zeros(T0+1,1);
Sf = zeros(T0+1,1);
cf = zeros(T0+1,1);
Q = zeros(T0+1,1);
dVk_f = zeros(T0+1,1);
dVS = zeros(T0+1,1);
gf = zeros(T0+1,1);
fval = zeros(T0+1,2);



Vf(1) = Vr(end);
kf(1) = kr(end);
cf(1) = cr(end);
invf(1) = ir(end);
Sf(1) = S_T0;

for i = 1:T0+1
    
    Q(i) = (1+popgr)^(T0-i+1);
    g = @(N) alpha0 + alpha1./(Sbar-alpha2./(alpha3+N)-Sf(i));
    gdS = @(N) alpha1*(alpha3+N).^2./((Sbar-Sf(i)).*(alpha3+N)-alpha2).^2;
    gdN = @(N) -alpha1.*alpha2./((Sbar-Sf(i)).*(alpha3+N)-alpha2).^2;
    
    if i == 1
        
        % The following function is only true when dV'/dS' = 0.
        funN = @(N) (1-g(N))*A-delta+A*kf(1)*gdN(N);
        Nf(i) = fzero(funN,43.5427);
        gf(i) = g(Nf(i));
        if gf(i) < alpha0
            error('N_T0 is too small given the guess of S_T0')
            
        end
        
        dVk_f(i) = cf(i)^(-gamma)*(1-A*kf(i)*gdN(Nf(i)));
        dVS(i) = -cf(i)^(-gamma)*A*kf(i)*gdS(Nf(i));
        
    else
        
        x0 = [kf(i-1),Nf(i-1)];
        
        options = optimset('Display','off','TolFun',1e-8,'TolX',1e-8);
        [xopt,fval(i,:)] = fsolve(@FOC_fossil_VF,x0,options);
        kf(i) = xopt(1);
        Nf(i) = xopt(2);
        
%         cf(i) = (beta*(1-A*kf(i-1)*gdNf(i-1))).^(-gamr)*cf(i-1);
        cf(i) = (beta*dVk_f(i-1)).^(-1/gamma);
        nf(i) = Nf(i-1)-Nf(i);
        invf(i) = kf(i-1)-(1-delta)*kf(i);
        Sf(i) = Sf(i-1)-Q(i)*A*kf(i);
        gf(i) = g(Nf(i));
        dVk_f(i) = cf(i)^(-gamma)*(A-g(Nf(i))*A+1-delta)+beta*dVS(i-1)*Q(i)*A;
        dVS(i) = -cf(i)^(-gamma)*A*kf(i)*gdS(Nf(i))+beta*dVS(i-1);
        
      
    end
    
    if Nf(i) < -1e-4
        error(['When i=', num2str(i), ',N=',num2str(Nf(i)),', S=',num2str(Sf(i)), ', k=',num2str(kf(i))],'Decrease Tbar or k_Tbar.');
        
    end
    if Sf(i) < 0
        error(['When i=', num2str(i), ',N=',num2str(Nf(i)),', S=',num2str(Sf(i)), ', k=',num2str(kf(i))]);
        
    end
end

disp(['k0 = ', num2str(kf(end)), ',N0 = ',num2str(Nf(end)),', S0 = ',num2str(Sf(end))]);

% Compare the results of the two models.
figure(2)
subplot(2,3,1), plot(kf(end:-1:1))
hold on
plot (kf0,'red')
title('k')
subplot(2,3,2), plot(Nf(end:-1:1))
hold on
plot (Nf0,'red')
title('N')
subplot(2,3,3), plot(Sf(end:-1:1))
hold on
plot (Sf0,'red')
title('S')
subplot(2,3,4), plot(invf(end:-1:1))
hold on
plot (if0,'red')
title('i')
subplot(2,3,5), plot(nf(end:-1:1))
hold on
plot (nf0,'red')
title('n')
subplot(2,3,6), plot(gf(end:-1:1))
hold on
plot (gf0,'red')
title('g')
toc
